Symbiodiniaceae community dynamics in Pocillopora colonies from Uva Island reef during 2014-2016
Ana Palacio
October 7, 2019
- 1. Symbiont to host (S/H) cell ratio and Durusdinium proportion
- 2. qPCR DATA CLEANING and quality control
- 3. Field/sample information
- 4. Community classification by genera presence/dominance
- 5a. qPCR based analysis (Symbiodiniaceae genera hosted)
- 6a. Changes in colony proportions based on detected genera (qPCR)
- 7a. Changes in total Symbiodiniaceae abundance (bleaching): qPCR
- 8a. Changes in Cladocopium abundance
- 9a. Changes in Durusdinium abundance
- 5b. DGGE based analysis (Symbiodiniaceae ITS2 type hosted)
- 6b. Changes in colony proportions based on ITS2 profiles (2014) and genus dominance
- 7b. Changes in total Symbiodiniaceae abundance (bleaching): qPCR + Initial ITS2
- 8b. Changes in Cladocopium abundance
- 9b. Changes in Durusdinium abundance
- 10. Supplmenmtary C/H D/H dynamics
This file analyzes the changes in S/H cell ratio (bleaching) in Pocillopora sp: * In 29 colonies sampled in 2014 (Before ENSO), 2015 (during 1st peak of ENSO) and 2016 (during 2nd peak of ENSO) * These samples were collected at Uva Island, Gulf of ChiriquÃ, Panama
# Libraries
library(tidyverse) # install.packages('tidyverse')
# library(devtools) # install.packages("devtools")
# devtools::install_github("jrcunning/steponeR")
library(steponeR)
library(scales)
library(gridExtra)
library(lme4)
library(lmerTest)
library(multcomp)
library(emmeans)
# Default ggplot settings
ggthe_bw<-theme(plot.background=element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)+
theme_bw()1. Symbiont to host (S/H) cell ratio and Durusdinium proportion
Get the raw data for all Pocillopora samples applying R.Cunning steponeR function
- Get list of plate files to read data
- Calculate the ratios
- Export results
Get colony and year information from sample labels
# Parse sample names, times, treatments, colonies, plates
sample.year <- rbind.fill(lapply(strsplit(as.character(Pdam$Sample.Name), split="_"),
function(X) data.frame(t(X))))
colnames(sample.year) <- c("Spp", "Colony","Year")
Pdam <- cbind(sample.year, Pdam[,-1])
# Keep unique sample ID to check reruns
Pdam$Sample<-paste(Pdam$Year,Pdam$Colony, sep="_")
Pdam$ID<-paste(Pdam$Sample,Pdam$File.Name, sep=".")Calculate total S/H ratio and D/C ratio and propD
Most of the colonies are dominated by one symbiont genus
2. qPCR DATA CLEANING and quality control
## 1. Check and remove NTC wells
ntc <- Pdam[which(Pdam$Colony=="NTC"), ]
Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(ntc), ])
## 2. Check and remove NoHSratio samples
NoHSratio <- Pdam[which(Pdam$tot.SH==0), ]
Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(NoHSratio), ])
## 3. Chose bw samples ran more than once
ReRunA <- Pdam[duplicated(Pdam$Sample),]
n_RunA <- data.frame(table(Pdam$Sample))
colnames(n_RunA)<-c("Sample","RanA")
Pdam<-join(Pdam, n_RunA, type = "left")
DuplicatesA <- Pdam[(Pdam$RanA>1),]
# Remove bad replicates
ToRem1<-read.csv("ToRemove10-28-16.csv") # 11/1/16
Pdam<-Pdam[!(Pdam$ID %in% ToRem1$ID),]
n_RunB <- data.frame(table(Pdam$Sample))
ReRunB <- Pdam[duplicated(Pdam$Sample),]
colnames(n_RunB)<-c("Sample","RanB")
Pdam<-join(Pdam, n_RunB, type = "left")
# List of dupplicated samples, should have 0 rows now
DuplicatesB <- Pdam[(Pdam$RanB>1),]
# 4. Check and remove samples with high SD or late amplification
StDe1.5 <- Pdam[which(Pdam$Pdam.CT.sd>1.5), ]
Pdam<-Pdam[!(Pdam$ID %in% StDe1.5$ID),]
# End of qPCR data cleaning3. Field/sample information
Get collection information and select only samples from Uva Island that were sampled in 2014, 2015 and 2016
# 1. Labeling and Factors
## Import treatment info from Field file
Info<-read.csv("SampleInfo.csv", header=TRUE)
## Merge SH ratios and sample info
Pdam<-join(Pdam, Info, by = "Sample",
type = "left", match = "all")
summary(Pdam)## Colony Year C.SH D.SH
## 256-1 : 3 14:142 Min. :0.000000 Min. :0.000000
## 257-1 : 3 15: 86 1st Qu.:0.000000 1st Qu.:0.000000
## 259-1 : 3 16: 56 Median :0.006687 Median :0.003965
## 260-1 : 3 Mean :0.070606 Mean :0.034897
## 262-1 : 3 3rd Qu.:0.036245 3rd Qu.:0.035000
## 264-1 : 3 Max. :2.046738 Max. :1.905453
## (Other):266
## Sample ID tot.SH
## Length:284 Length:284 Min. :0.001125
## Class :character Class :character 1st Qu.:0.018868
## Mode :character Mode :character Median :0.035563
## Mean :0.105503
## 3rd Qu.:0.078992
## Max. :2.046738
##
## logTot.SH logC.SH logD.SH D.Prp
## Min. :-2.9487 Min. :-5.5123 Min. :-5.726 Min. :0.0000
## 1st Qu.:-1.7243 1st Qu.:-2.4703 1st Qu.:-1.862 1st Qu.:0.0000
## Median :-1.4490 Median :-1.7240 Median :-1.544 Median :0.2061
## Mean :-1.3883 Mean :-1.9408 Mean :-1.871 Mean :0.4674
## 3rd Qu.:-1.1024 3rd Qu.:-1.2207 3rd Qu.:-1.250 3rd Qu.:1.0000
## Max. : 0.3111 Max. : 0.3111 Max. : 0.280 Max. :1.0000
## NA's :77 NA's :114
## RanA RanB Region Location
## Min. :1.000 Min. :1 Gal_Darwin : 58 Canal de afuera : 15
## 1st Qu.:1.000 1st Qu.:1 Gal_Floreana: 12 Contadora Island: 36
## Median :1.000 Median :1 Pan_Chiriqui:158 Darwin : 58
## Mean :1.092 Mean :1 Pan_Panama : 56 Floreana : 12
## 3rd Qu.:1.000 3rd Qu.:1 Saboga Island : 20
## Max. :3.000 Max. :1 Secas Island : 16
## Uva Island :127
## Reef Depht Spp
## Coral Collection Point : 40 Min. : 3.00 P.dam:164
## By the Pier : 34 1st Qu.:10.00 P.eff: 15
## Darwin Reef, next to BEAMS: 33 Median :20.00 P.ele: 61
## Shallow Millepora : 24 Mean :23.24 P.eyd: 14
## M.intrincata Shallow : 17 3rd Qu.:36.00 P.mea: 12
## Ridley Rock Shallow : 16 Max. :68.00 P.ver: 18
## (Other) :120
## Year.2 Number Tag Collection.day
## Min. :2014 Min. :1 no tag : 9 : 70
## 1st Qu.:2014 1st Qu.:1 270-1 : 4 08/19/14: 21
## Median :2014 Median :1 336 : 4 08/18/15: 17
## Mean :2015 Mean :1 339-1 : 4 08/27/15: 17
## 3rd Qu.:2015 3rd Qu.:1 340-1 : 4 04/18/16: 16
## Max. :2016 Max. :1 345 : 4 04/16/16: 14
## NA's :78 (Other):255 (Other) :129
## Notes Depth.2
## :274 mid : 86
## 1st Rock at 60ft, not the main rock: 3 Shallow:198
## No tag : 1
## No Tag : 3
## Not 259 : 1
## P.cap : 1
## Sampled 16 and 17 : 1# 2. Uva Data
# Write how many years a colony was sampled
Years <- data.frame(table(Pdam$Colony))
colnames(Years)<-c("Colony","Times")
Pdam<-join(Pdam, Years, type = "left")
# Select data only from Uva Island
data.Uva<-Pdam[which(Pdam$Location=="Uva Island"), ]
# Select colonies that were sampled three times (2014, 2015 and 2016)
data.Uva<-data.Uva[which(data.Uva$Times>2),]
# Remove "Ross" samples (colony in deep location (~50feet), not Uva reef)
# and 259 as well as 264 that were mislabeled in at least one of the sampling points
data.Uva<-filter(data.Uva, Colony!="Ross")
data.Uva<-filter(data.Uva, Colony!="259-1")
data.Uva<-filter(data.Uva, Colony!="264-1")
data.Uva<-data.Uva[,c("Colony","Year.2","C.SH","D.SH","tot.SH",
"logTot.SH","D.Prp", "logC.SH", "logD.SH","Year","Sample")]
data.Uva$Year<-factor(data.Uva$Year, levels=c("14", "15", "16"))
data.Uva$Colony<-factor(data.Uva$Colony)
summary(data.Uva)## Colony Year.2 C.SH D.SH
## 256-1 : 3 Min. :2014 Min. :0.000000 Min. :0.00000
## 257-1 : 3 1st Qu.:2014 1st Qu.:0.000000 1st Qu.:0.00000
## 260-1 : 3 Median :2015 Median :0.003005 Median :0.01442
## 262-1 : 3 Mean :2015 Mean :0.019174 Mean :0.02071
## 278-1 : 3 3rd Qu.:2016 3rd Qu.:0.011841 3rd Qu.:0.03001
## 284-1 : 3 Max. :2016 Max. :0.690021 Max. :0.18875
## (Other):69
## tot.SH logTot.SH D.Prp logC.SH
## Min. :0.00146 Min. :-2.8356 Min. :0.0000 Min. :-5.0627
## 1st Qu.:0.01294 1st Qu.:-1.8888 1st Qu.:0.0000 1st Qu.:-2.5001
## Median :0.02266 Median :-1.6448 Median :0.8620 Median :-2.0971
## Mean :0.03988 Mean :-1.6539 Mean :0.5411 Mean :-2.1808
## 3rd Qu.:0.03579 3rd Qu.:-1.4463 3rd Qu.:1.0000 3rd Qu.:-1.7015
## Max. :0.69002 Max. :-0.1611 Max. :1.0000 Max. :-0.1611
## NA's :32
## logD.SH Year Sample
## Min. :-5.7262 14:29 Length:87
## 1st Qu.:-1.8381 15:29 Class :character
## Median :-1.6245 16:29 Mode :character
## Mean :-2.0059
## 3rd Qu.:-1.4640
## Max. :-0.7241
## NA's :294. Community classification by genera presence/dominance
Symbiodiniaceae genera detected in the samples (Only-C, Only-D, Mixed communities)
# Add a column for the genera detected using qPCR
data.Uva$Community<-NULL
data.Uva$Community[which(data.Uva$D.Prp>=0)] <- "Mixed"
data.Uva$Community[which(data.Uva$D.Prp==0)] <- "Cladocopium-only"
data.Uva$Community[which(data.Uva$D.Prp==1)] <- "Durusdinium-only"
data.Uva$Community<-factor(as.character(data.Uva$Community),
levels=c("Cladocopium-only", "Mixed", "Durusdinium-only"))
# Summary genera detected
Summary_Community<-data.Uva %>%
group_by(Year.2, Community) %>%
dplyr::summarise(total.count=n())
Summary_CommunityDominant Symbiodiniaceae genus in each sample (C-dominated, D-dominated)
In general, Pocillipora colonies are dominated by one Symbiodiniaceae type that accounts for more than 90 of the symbiont community. However, sometimes intermediate proportions are common when the symbiont community is transitioning from being dominated by one algal type to other (e.g. during heat stress).
For this data set, 2 of the colonies sampled in 2014 had relativelly even proportions of Cladocopium_b and Durusdinim glynnii (Colony 267_2014, D proportion = 0.47; Colony 269_2014, D proportion = 0.50). This made difficult their classification into one domminant type for statistical testing and further plotting. Since these 2 colonies were later dominated by Durusdinim glynnii in 2015, we decided to include them in the D-dominated classification, eventhoug in 2014 there were not such a dominance (See figure 2).
# Add a column for the dominant genus
data.Uva$DominantCommunity<-NULL
#data.Uva$Community[which(data.Uva$D.Prp>0.1 & data.Uva$D.Prp<0.50)] <- "CD"
#data.Uva$Community[which(data.Uva$D.Prp>0.5 & data.Uva$D.Prp<0.9)] <- "DC"
data.Uva$DominantCommunity[which(data.Uva$D.Prp>=0.45)] <- "Durusdinium-dominated"
data.Uva$DominantCommunity[which(data.Uva$D.Prp<0.45)] <- "Cladocopium-dominated"
data.Uva$DominantCommunity<-factor(as.character(data.Uva$DominantCommunity),
levels=c("Cladocopium-dominated","Durusdinium-dominated"))
# Summary dominant genera
Summary_Community2<-data.Uva %>%
group_by(Year.2, DominantCommunity) %>%
dplyr::summarise(total.count=n())
Summary_Community2Changes in the number of colonies dominated by each genus
To test the statistical significance of changes in the number of colonies dominated by Durusdinium versus Cladocopium we considered inclussion and exclussion of these 2 initially even colonies, and to place them in the D-dominated or C-dominated categories in 2014.
- Fisher test including initially (Aug 2014) even colonies. Both colonies placed in the D-dominated category.
Input =("
Year C D
2014 15 14
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.1064
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8336155 9.7658946
## sample estimates:
## odds ratio
## 2.760908- Fisher test including initially (Aug 2014) even colonies. One colony was placed in the C-dominated category (D proportion=0.47), while the second one was placed in the D-dominated category (D proportion = 0.50)
Input =("
Year C D
2014 16 13
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.06099
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.9556883 11.2396870
## sample estimates:
## odds ratio
## 3.162767- Fisher test removing initially (Aug 2014) even colonies. N=27 colonies instead 29
Input =("
Year C D
2014 15 12
2016 8 19
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.09778
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.848568 10.659493
## sample estimates:
## odds ratio
## 2.906792- Fisher test including initially (Aug 2014) even colonies. Both colonies placed in the C-dominated category (D proportion=0.47, colony 267; and 0.50, colony 269)
Input =("
Year C D
2014 17 12
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.03295
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.094473 13.006776
## sample estimates:
## odds ratio
## 3.6297775a. qPCR based analysis (Symbiodiniaceae genera hosted)
Set initial (2014) and final (2016) communities for each sample
# Subset Aug 2014, Aug 2015 and April 2016 samples
D2014.data<-subset(data.Uva, Year.2==2014)
D2014.data<-D2014.data[,c("Colony","D.Prp","Community","DominantCommunity")]
colnames(D2014.data)<-(c("Colony","D14.D.Prp","D14Community", "D14Dominant"))
data.Uva<-left_join(data.Uva, D2014.data, by=c("Colony")) # Add initial community to all samples
# D2015.data<-subset(data.Uva, Year.2==2015)
# D2015.data<-D2015.data[, c("Colony","D.Prp","Community")]
# colnames(D2015.data)<-(c("Colony","D15.D.Prp","D15Community"))
# # data.Uva<-left_join(data.Uva, D2015.data, by=c("Colony"))
D2016.data<-subset(data.Uva, Year.2==2016) # D1a under control temperature
D2016.data<-D2016.data[, c("Colony","D.Prp","Community","DominantCommunity")]
colnames(D2016.data)<-(c("Colony","D16.D.Prp","D16Community", "D16Dominant"))
data.Uva<-left_join(data.Uva, D2016.data, by=c("Colony")) # Add final community to all samples6a. Changes in colony proportions based on detected genera (qPCR)
- Based on detected Symbiodiniacea genera in 2014 (qPCR)
Histo_Dprop <- ggplot(data.Uva, aes (D.Prp , fill=(D14Community))) + ggthe_bw +
geom_histogram(binwidth = 0.05, boundary=0) +
facet_grid(~Year.2, labeller = labeller(Year.2=c(
`2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~(D/H)/(S/H)),
breaks=(seq(0,1,by=0.1))) +
scale_y_continuous((name= "Number of colonies"),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
Histo_Dprop +coord_flip()
Traject_Proportion_Facet <- ggplot (data.Uva) +
ggthe_bw +
geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+
facet_wrap(~D14Community, labeller = labeller(
D14Community=c(`Cladocopium-only` = "Cladocopium-only (2014)",
`Mixed` = "Mixed (2014)",
`Durusdinium-only` = "Durusdinium-only (2014)"))) +
geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3,
height = 0.05, width = 0.15) +
#labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
theme(legend.position="none")
Traject_Proportion_Facet
Figure 2: qPCR
Histo2014_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(D14Community))) +
ggthe_bw + coord_flip()+ ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2014_Dprop
Histo2016_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(D14Community))) +
ggthe_bw + coord_flip()+ ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2016_Dprop
Traject_Proportion<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , colour=(D14Community), group=(Colony))) + ggthe_bw + ggtitle("b.") +
geom_line(linetype=2) +
geom_jitter(height = 0.01, width = 0.25, alpha=0.4) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.01,0.01)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Traject_Proportion
Figure_2<-grid.arrange(Histo2014_Dprop,Traject_Proportion, Histo2016_Dprop,
nrow=1, widths=c(3/9, 3/9, 3/9))
# ggsave(file="Outputs/Figure_2qPCR.svg", plot=Figure_2, width=7, height=4)7a. Changes in total Symbiodiniaceae abundance (bleaching): qPCR
- Only considering initial (2014) community
SH_2014Co<-lmerTest::lmer(logTot.SH ~ Year * D14Community + (1|Colony), data=data.Uva)
summary(SH_2014Co)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 83.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6538 -0.5202 -0.1347 0.3700 2.9430
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.03036 0.1742
## Residual 0.10764 0.3281
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.6402 0.1238 71.1153 -13.246
## Year15 -0.4132 0.1547 52.0000 -2.672
## Year16 -0.5536 0.1547 52.0000 -3.579
## D14CommunityMixed -0.0137 0.1638 71.1153 -0.084
## D14CommunityDurusdinium-only 0.1232 0.1805 71.1153 0.682
## Year15:D14CommunityMixed 0.7731 0.2046 52.0000 3.778
## Year16:D14CommunityMixed 0.6960 0.2046 52.0000 3.402
## Year15:D14CommunityDurusdinium-only 0.4584 0.2255 52.0000 2.033
## Year16:D14CommunityDurusdinium-only 0.3869 0.2255 52.0000 1.716
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.010053 *
## Year16 0.000756 ***
## D14CommunityMixed 0.933599
## D14CommunityDurusdinium-only 0.497283
## Year15:D14CommunityMixed 0.000408 ***
## Year16:D14CommunityMixed 0.001294 **
## Year15:D14CommunityDurusdinium-only 0.047152 *
## Year16:D14CommunityDurusdinium-only 0.092125 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 D14CmM D14CD- Y15:D14CM Y16:D14CM
## Year15 -0.624
## Year16 -0.624 0.500
## D14CmmntyMx -0.756 0.472 0.472
## D14CmmntyD- -0.686 0.428 0.428 0.519
## Yr15:D14CmM 0.472 -0.756 -0.378 -0.624 -0.324
## Yr16:D14CmM 0.472 -0.378 -0.756 -0.624 -0.324 0.500
## Yr15:D14CD- 0.428 -0.686 -0.343 -0.324 -0.624 0.519 0.259
## Yr16:D14CD- 0.428 -0.343 -0.686 -0.324 -0.624 0.259 0.519
## Y15:D14CD
## Year15
## Year16
## D14CmmntyMx
## D14CmmntyD-
## Yr15:D14CmM
## Yr16:D14CmM
## Yr15:D14CD-
## Yr16:D14CD- 0.500 #step(SH_2014Co)
SH.emmc<-emmeans(SH_2014Co, ~Year*D14Community)
SH_groups<-multcomp::cld(SH.emmc)
SH_groups plot(emmeans(SH_2014Co, ~Year|~D14Community), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~D14Community)
SH_2014Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(D14Community))) + ggthe_bw +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position=position_dodge(width=0.95))+
stat_summary(fun.data = "mean_cl_boot",geom = "point",
position=position_dodge(width=0.95))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Symbionts detectec in 2014",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.5,0.1),
strip.background =element_rect(fill="white")) +
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_discrete(name="Year")
SH_2014Com
- Only considering current community
SH_Community<-lmerTest::lmer(logTot.SH ~ Year * Community + (1|Colony), data=data.Uva)
summary(SH_Community)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 91.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3973 -0.5171 -0.0987 0.3905 4.0587
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02404 0.1550
## Residual 0.12526 0.3539
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.60009 0.12804 76.71005 -12.497
## Year15 -0.20778 0.15540 49.73820 -1.337
## Year16 -0.71036 0.17979 47.25496 -3.951
## CommunityMixed -0.09453 0.16829 77.96039 -0.562
## CommunityDurusdinium-only 0.09882 0.18645 77.17623 0.530
## Year15:CommunityMixed 0.58170 0.24003 54.72488 2.423
## Year16:CommunityMixed 0.81742 0.24650 54.96059 3.316
## Year15:CommunityDurusdinium-only 0.27313 0.22958 47.93524 1.190
## Year16:CommunityDurusdinium-only 0.60594 0.24059 49.84009 2.519
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.187280
## Year16 0.000258 ***
## CommunityMixed 0.575911
## CommunityDurusdinium-only 0.597640
## Year15:CommunityMixed 0.018708 *
## Year16:CommunityMixed 0.001622 **
## Year15:CommunityDurusdinium-only 0.240035
## Year16:CommunityDurusdinium-only 0.015039 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmntM CmmnD- Y15:CM Y16:CM Y15:CD
## Year15 -0.733
## Year16 -0.602 0.497
## CommuntyMxd -0.754 0.576 0.457
## CmmntyDrsd- -0.686 0.504 0.414 0.524
## Yr15:CmmntM 0.472 -0.658 -0.322 -0.628 -0.326
## Yr16:CmmntM 0.459 -0.365 -0.740 -0.608 -0.318 0.412
## Yr15:CmmnD- 0.496 -0.677 -0.337 -0.383 -0.717 0.440 0.243
## Yr16:CmmnD- 0.451 -0.369 -0.747 -0.326 -0.682 0.240 0.543 0.541 #step(SH_Community)
SH.emmc<-emmeans(SH_Community, ~Year*Community)
SH_groups<-multcomp::cld(SH.emmc)
#write.csv(SH_groups,"SH_groupsB.csv")
plot(emmeans(SH_Community, ~Year|~Community), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~Community)
SH_Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position=position_dodge(width=0.95))+
stat_summary(fun.data = "mean_cl_boot",geom = "point",
position=position_dodge(width=0.95))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Symbionts detectec in 2014",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.5,0.1),
strip.background =element_rect(fill="white")) +
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_discrete(name="Year")
SH_Com
anova(SH_2014Co, SH_Community)- Plot based on initial and current symbiont community
Figure 3: qPCR
SH_Communities<-lmerTest::lmer(logTot.SH ~ Year * D14Community * Community + (1|Colony), data=data.Uva)
summary(SH_Communities)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community * Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 81.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73407 -0.48463 -0.06579 0.35339 2.91935
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02179 0.1476
## Residual 0.11026 0.3320
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.6402 0.1211 70.1550 -13.541
## Year15 -0.4132 0.1565 47.2363 -2.640
## Year16 -0.7093 0.1688 50.6884 -4.202
## D14CommunityMixed -0.7146 0.3265 73.0936 -2.189
## D14CommunityDurusdinium-only -0.6296 0.3928 72.9975 -1.603
## CommunityMixed 0.7009 0.2845 70.5547 2.464
## CommunityDurusdinium-only 0.7527 0.3508 70.9827 2.146
## Year15:D14CommunityMixed 1.5019 0.3806 64.8657 3.946
## Year16:D14CommunityMixed 0.8259 0.2396 53.9251 3.447
## Year15:D14CommunityDurusdinium-only 1.2674 0.5110 61.4046 2.480
## Year16:D14CommunityDurusdinium-only 0.5426 0.2368 48.9752 2.292
## Year15:CommunityMixed -0.7379 0.3653 70.6556 -2.020
## Year15:CommunityDurusdinium-only -0.8090 0.4572 64.7392 -1.770
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.011206 *
## Year16 0.000107 ***
## D14CommunityMixed 0.031805 *
## D14CommunityDurusdinium-only 0.113266
## CommunityMixed 0.016178 *
## CommunityDurusdinium-only 0.035332 *
## Year15:D14CommunityMixed 0.000198 ***
## Year16:D14CommunityMixed 0.001108 **
## Year15:D14CommunityDurusdinium-only 0.015876 *
## Year16:D14CommunityDurusdinium-only 0.026258 *
## Year15:CommunityMixed 0.047162 *
## Year15:CommunityDurusdinium-only 0.081500 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 14 columns / coefficients #step(SH_Communities)
SH.emmc<-emmeans(SH_Communities, ~Year*D14Community*Community)
SH_groups<-cld(SH.emmc)
SH_groups<-as.data.frame(SH_groups[complete.cases(SH_groups),])
SH_groups<-SH_groups[order(SH_groups$Year,SH_groups$D14Community, SH_groups$Community),]
SH_groups #write.csv(SH_groups,"SH_Communities_multicomp.csv")
SH_Comties <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw+
facet_grid(~D14Community, labeller = labeller(D14Community=c(
`Cladocopium-only` = "Cladocopium-only (2014)",
`Mixed` = "Mixed (2014)",
`Durusdinium-only` = "Durusdinium-only (2014)"))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.2),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(14, 15, 16),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
SH_Comties
#ggsave(file="Figure_3.svg", plot=SH_Comties, width=6, height=4)8a. Changes in Cladocopium abundance
CH_Com <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(Community))) + ggthe_bw+ ggtitle("a.")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.8),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.7),
name="\n Cladocopium to host cell ratio (C/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" ", " ", " "))
#CH_Com
CH_Com2<-CH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
`Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
`Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))
CH<-lmer(logC.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
summary(CH)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * Community * D14Dominant + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 113.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70311 -0.36865 -0.09556 0.50582 1.89250
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.1182 0.3438
## Residual 0.4032 0.6350
## Number of obs: 55, groups: Colony, 21
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -1.54068 0.23810 44.78323
## Year15 -0.25625 0.28031 26.58185
## Year16 -0.73614 0.32339 24.71137
## CommunityMixed -0.32887 0.37022 44.99572
## D14DominantDurusdinium-dominated -1.09559 0.41353 43.96389
## Year15:CommunityMixed 0.59388 0.61658 32.70059
## Year16:CommunityMixed -0.03386 0.50241 28.81910
## Year15:D14DominantDurusdinium-dominated -0.54917 0.68606 30.91184
## Year16:D14DominantDurusdinium-dominated 0.88798 0.66184 31.43058
## t value Pr(>|t|)
## (Intercept) -6.471 6.36e-08 ***
## Year15 -0.914 0.3688
## Year16 -2.276 0.0318 *
## CommunityMixed -0.888 0.3791
## D14DominantDurusdinium-dominated -2.649 0.0112 *
## Year15:CommunityMixed 0.963 0.3425
## Year16:CommunityMixed -0.067 0.9467
## Year15:D14DominantDurusdinium-dominated -0.800 0.4296
## Year16:D14DominantDurusdinium-dominated 1.342 0.1893
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmntM D14DD- Y15:CM Y16:CM Y15:D1
## Year15 -0.715
## Year16 -0.579 0.492
## CommuntyMxd -0.622 0.502 0.369
## D14DmnntDr- -0.019 -0.038 0.003 -0.537
## Yr15:CmmntM 0.318 -0.474 -0.224 -0.512 0.275
## Yr16:CmmntM 0.400 -0.331 -0.659 -0.643 0.345 0.357
## Yr15:D14DD- 0.006 0.017 0.000 0.255 -0.468 -0.705 -0.185
## Yr16:D14DD- -0.021 0.011 0.011 0.308 -0.509 -0.162 -0.437 0.306
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients #step(CH)
CH.emmc<-emmeans(CH, ~Year*Community* D14Dominant)
CH_groups<-cld(CH.emmc)
CH_groups<-as.data.frame(CH_groups[complete.cases(CH_groups),])
CH_groups<-CH_groups[order(CH_groups$Year,CH_groups$D14Dominant, CH_groups$Community),]
CH_groups #write.csv(CH_groups,"CH_groups.csv")9a. Changes in Durusdinium abundance
DH_Com <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(Community))) +
ggthe_bw+ ggtitle("b.")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.2),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-6, -0.7),
name="\n Durusdinium to host cell ratio (D/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
# DH_Com
DH_Com2<-DH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
`Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
`Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))
DH_LM<-lmer(logD.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
summary(DH_LM)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * Community * D14Dominant + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 66.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0953 -0.5756 0.0079 0.3844 3.6310
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.0000 0.0000
## Residual 0.1687 0.4107
## Number of obs: 58, groups: Colony, 22
##
## Fixed effects:
## Estimate
## (Intercept) -4.86573
## Year15 0.07049
## Year16 3.21459
## CommunityDurusdinium-only 0.27806
## D14DominantDurusdinium-dominated 3.15614
## Year15:CommunityDurusdinium-only -0.35506
## Year16:CommunityDurusdinium-only -0.22094
## Year15:D14DominantDurusdinium-dominated 0.36688
## Year16:D14DominantDurusdinium-dominated -3.09697
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated -0.08557
## Std. Error
## (Intercept) 0.16766
## Year15 0.33532
## Year16 0.23711
## CommunityDurusdinium-only 0.50992
## D14DominantDurusdinium-dominated 0.23711
## Year15:CommunityDurusdinium-only 0.32897
## Year16:CommunityDurusdinium-only 0.38416
## Year15:D14DominantDurusdinium-dominated 0.42745
## Year16:D14DominantDurusdinium-dominated 0.41068
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 0.45916
## df
## (Intercept) 48.00000
## Year15 48.00000
## Year16 48.00000
## CommunityDurusdinium-only 48.00000
## D14DominantDurusdinium-dominated 48.00000
## Year15:CommunityDurusdinium-only 48.00000
## Year16:CommunityDurusdinium-only 48.00000
## Year15:D14DominantDurusdinium-dominated 48.00000
## Year16:D14DominantDurusdinium-dominated 48.00000
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 48.00000
## t value
## (Intercept) -29.021
## Year15 0.210
## Year16 13.557
## CommunityDurusdinium-only 0.545
## D14DominantDurusdinium-dominated 13.311
## Year15:CommunityDurusdinium-only -1.079
## Year16:CommunityDurusdinium-only -0.575
## Year15:D14DominantDurusdinium-dominated 0.858
## Year16:D14DominantDurusdinium-dominated -7.541
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated -0.186
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.834
## Year16 < 2e-16 ***
## CommunityDurusdinium-only 0.588
## D14DominantDurusdinium-dominated < 2e-16 ***
## Year15:CommunityDurusdinium-only 0.286
## Year16:CommunityDurusdinium-only 0.568
## Year15:D14DominantDurusdinium-dominated 0.395
## Year16:D14DominantDurusdinium-dominated 1.09e-09 ***
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 0.853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmnD- D14DD- Y15:CD Y16:CD Y15:D1 Y16:D1
## Year15 -0.500
## Year16 -0.707 0.354
## CmmntyDrsd- 0.000 0.000 -0.232
## D14DmnntDr- -0.707 0.354 0.500 -0.232
## Yr15:CmmnD- 0.000 0.000 0.000 -0.293 0.360
## Yr16:CmmnD- 0.000 0.000 0.000 -0.753 0.309 0.389
## Yr15:D14DD- 0.392 -0.784 -0.277 0.129 -0.555 -0.500 -0.171
## Yr16:D14DD- 0.408 -0.204 -0.577 0.671 -0.577 -0.208 -0.713 0.320
## CmD-:D14DD- 0.000 0.000 0.258 -0.900 0.000 0.000 0.558 0.000 -0.596
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## convergence code: 0
## boundary (singular) fit: see ?isSingular #step(DH_LM)
DH.emmc<-emmeans(DH_LM, ~Year*Community* D14Dominant)
DH_groups<-cld(DH.emmc)
DH_groups<-as.data.frame(DH_groups[complete.cases(DH_groups),])
DH_groups<-DH_groups[order(DH_groups$Year,DH_groups$D14Dominant, DH_groups$Community),]
DH_groups #write.csv(DH_groups,"DH_groups.csv")Figure 4: qPCR
Figure_4<-grid.arrange(CH_Com2, DH_Com2,nrow=2)
# ggsave(file="Outputs/Figure_4.svg", plot=Figure_4, width=7, height=4)5b. DGGE based analysis (Symbiodiniaceae ITS2 type hosted)
Community classification by ITS2 dominant profile in 2014
# Add initial (2014) ITS2 profiles
ITS2_C_a<-data.frame("Colony"=c("256-1", "257-1", "260-1", "262-1", "330-1", "340-1A"),
"ITS2_14"="C_a")
# ITS2_C_b<-data.frame("Colony"=c("265-1", "266-1", "267", "268-1", "269-1", "270-1A", "271-1", "285-1", "287-1", "336-1A", "339-1A"), "ITS2_14"="C_b") # 267 and 269 placed under D1
ITS2_C_b<-data.frame("Colony"=c("265-1", "266-1", "268-1", "270-1A", "271-1", "285-1",
"287-1", "336-1A", "339-1A"), "ITS2_14"="C_b")
ITS2_C<-rbind(ITS2_C_a, ITS2_C_b)
ITS2_D<-anti_join(D2014.data, ITS2_C, by="Colony")
ITS2_D$ITS2_14<-"D1"
ITS2_D<-ITS2_D[, c("Colony", "ITS2_14")]
ITS2<-rbind(ITS2_C, ITS2_D)
data.Uva<-left_join(data.Uva, ITS2, by=c("Colony")) # Add initial ITS2 profiles to all samples6b. Changes in colony proportions based on ITS2 profiles (2014) and genus dominance
- Based on ITS2 detected symbionts in 2014
Histo_Dprop_ITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(ITS2_14)), shape=(D14Community)) + ggthe_bw +
geom_histogram(binwidth = 0.05, boundary=0) +
facet_grid(~Year.2, labeller = labeller(Year.2=c(
`2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("C_a", "C_b", "D1")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~(D/H)/(S/H)),
breaks=(seq(0,1,by=0.1))) +
scale_y_continuous((name= "Number of colonies"),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
Histo_Dprop_ITS2 +coord_flip()
Traject_Proportion_FacetITS2 <- ggplot (data.Uva) +
ggthe_bw +
geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+
facet_wrap(~ITS2_14) +
geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3,
height = 0.05, width = 0.15) +
#labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
theme(legend.position="none")
Traject_Proportion_FacetITS2
Figure 2: ITS2
Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(ITS2_14))) +
ggthe_bw + coord_flip()+ ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("C_a", "C_b", "D1")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2014_DpropITS2
Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(ITS2_14))) +
ggthe_bw + coord_flip()+ ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("C_a", "C_b", "D1")) +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2016_DpropITS2
Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , colour=(ITS2_14), group=(Colony)))+
ggthe_bw + ggtitle("b.") +
geom_line(linetype=2) +
geom_jitter(height = 0.01, width = 0.25, alpha=0.4) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("C_a", "C_b", "D1")) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.01,0.01)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Traject_ProportionITS2
Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(3/9, 3/9, 3/9))
# ggsave(file="Outputs/Figure_ITS2.svg", plot=Figure_2_ITS2, width=7, height=4)7b. Changes in total Symbiodiniaceae abundance (bleaching): qPCR + Initial ITS2
Only considering initial (2014) community
SH_2014CoI<-lmerTest::lmer(logTot.SH ~ Year * ITS2_14 + (1|Colony), data=data.Uva)
summary(SH_2014CoI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * ITS2_14 + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 77.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3417 -0.5054 -0.0500 0.3438 3.4861
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.0190 0.1378
## Residual 0.1058 0.3253
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.66014 0.14424 74.54698 -11.510 < 2e-16 ***
## Year15 -0.56741 0.18782 52.00000 -3.021 0.003902 **
## Year16 -0.77804 0.18782 52.00000 -4.142 0.000127 ***
## ITS2_14C_b -0.01746 0.18621 74.54698 -0.094 0.925547
## ITS2_14D1 0.11106 0.17240 74.54698 0.644 0.521417
## Year15:ITS2_14C_b 0.79408 0.24248 52.00000 3.275 0.001885 **
## Year16:ITS2_14C_b 0.88586 0.24248 52.00000 3.653 0.000602 ***
## Year15:ITS2_14D1 0.73349 0.22449 52.00000 3.267 0.001926 **
## Year16:ITS2_14D1 0.71310 0.22449 52.00000 3.177 0.002508 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 ITS2_14C ITS2_14D Y15:ITS2_14C
## Year15 -0.651
## Year16 -0.651 0.500
## ITS2_14C_b -0.775 0.504 0.504
## ITS2_14D1 -0.837 0.545 0.545 0.648
## Y15:ITS2_14C 0.504 -0.775 -0.387 -0.651 -0.422
## Y16:ITS2_14C 0.504 -0.387 -0.775 -0.651 -0.422 0.500
## Y15:ITS2_14D 0.545 -0.837 -0.418 -0.422 -0.651 0.648
## Y16:ITS2_14D 0.545 -0.418 -0.837 -0.422 -0.651 0.324
## Y16:ITS2_14C Y15:ITS2_14D
## Year15
## Year16
## ITS2_14C_b
## ITS2_14D1
## Y15:ITS2_14C
## Y16:ITS2_14C
## Y15:ITS2_14D 0.324
## Y16:ITS2_14D 0.648 0.500 #step(SH_2014CoI)
SH_I.emmc<-emmeans(SH_2014CoI, ~Year*ITS2_14)
SH_I_groups<-multcomp::cld(SH_I.emmc)
SH_I_groups SH_I_groups<-as.data.frame(SH_I_groups[complete.cases(SH_I_groups),])
SH_I_groups<-SH_I_groups[order(SH_I_groups$Year,SH_I_groups$ITS2_14),]
SH_I_groups plot(emmeans(SH_2014CoI, ~Year|~ITS2_14), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~ITS2_14)
SH_2014ComI <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(ITS2_14))) + ggthe_bw +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position=position_dodge(width=0.95))+
stat_summary(fun.data = "mean_cl_boot",geom = "point",
position=position_dodge(width=0.95))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "ITS2 detected in 2014",
labels=c("C_a", "C_b", "D1"))+
theme(legend.position=c(0.5,0.1),
strip.background =element_rect(fill="white")) +
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_discrete(name="Year")
SH_2014ComI
Only considering current community ADD CURRENT ITS2!
# SH_Community<-lmerTest::lmer(logTot.SH ~ Year * Community + (1|Colony), data=data.Uva)
# summary(SH_Community)
# #step(SH_Community)
#
# SH.emmc<-emmeans(SH_Community, ~Year*Community)
# SH_groups<-multcomp::cld(SH.emmc)
# #write.csv(SH_groups,"SH_groupsB.csv")
#
# plot(emmeans(SH_Community, ~Year|~Community), comparisons = TRUE) +
# coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~Community)
#
# SH_Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw +
# stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
# position=position_dodge(width=0.95))+
# stat_summary(fun.data = "mean_cl_boot",geom = "point",
# position=position_dodge(width=0.95))+
# stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
#
# scale_colour_manual(values=c("blue","purple", "red"),
# guide = guide_legend(title.position = "top", nrow = 1),
# name = "Symbionts detectec in 2014",
# labels=c("Cladocopium", "Mixed", "Durusdinium"))+
# theme(legend.position=c(0.5,0.1),
# strip.background =element_rect(fill="white")) +
# scale_y_continuous(limits = c(-3, -0.8),
# name="\n Total symbiont to host cell ratio (S/H)",
# labels = (math_format(10^.x)),
# expand = c(0.0, 0)) +
# scale_x_discrete(name="Year")
# SH_Com- Plot based on initial and current symbiont community
Figure 3: ITS2
SH_Communities_I<-lmerTest::lmer(logTot.SH ~ Year * ITS2_14 * Community + (1|Colony), data=data.Uva)
summary(SH_Communities_I)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * ITS2_14 * Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 80.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2665 -0.4690 -0.1486 0.3796 3.3861
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02122 0.1457
## Residual 0.11122 0.3335
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.66014 0.14857 67.55796 -11.174
## Year15 -0.56741 0.19255 44.59436 -2.947
## Year16 -0.77804 0.19255 44.59436 -4.041
## ITS2_14C_b 0.10119 0.25532 70.67197 0.396
## ITS2_14D1 0.20013 0.51863 70.63905 0.386
## CommunityMixed -0.15133 0.47929 70.99075 -0.316
## CommunityDurusdinium-only -0.04237 0.50922 68.89495 -0.083
## Year15:ITS2_14C_b 0.67029 0.30501 51.27088 2.198
## Year16:ITS2_14C_b 0.44786 0.44197 54.72157 1.013
## Year15:ITS2_14D1 0.72446 0.47149 54.95813 1.537
## Year16:ITS2_14D1 0.24593 0.56263 57.15914 0.437
## Year15:CommunityMixed 0.20111 0.37100 59.52765 0.542
## Year16:CommunityMixed 0.53106 0.44747 59.20686 1.187
## Year15:CommunityDurusdinium-only -0.08959 0.45989 56.55615 -0.195
## Year16:CommunityDurusdinium-only 0.42021 0.55196 58.74860 0.761
## ITS2_14C_b:CommunityMixed -0.02664 0.39946 69.55482 -0.067
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.005090 **
## Year16 0.000208 ***
## ITS2_14C_b 0.693066
## ITS2_14D1 0.700746
## CommunityMixed 0.753132
## CommunityDurusdinium-only 0.933928
## Year15:ITS2_14C_b 0.032520 *
## Year16:ITS2_14C_b 0.315370
## Year15:ITS2_14D1 0.130144
## Year16:ITS2_14D1 0.663686
## Year15:CommunityMixed 0.589788
## Year16:CommunityMixed 0.240040
## Year15:CommunityDurusdinium-only 0.846244
## Year16:CommunityDurusdinium-only 0.449526
## ITS2_14C_b:CommunityMixed 0.947010
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 11 columns / coefficients #step(SH_Communities_I)
SH_CI.emmc<-emmeans(SH_Communities_I, ~Year*ITS2_14*Community)
SH_CI_groups<-cld(SH_CI.emmc)
SH_CI_groups<-as.data.frame(SH_CI_groups[complete.cases(SH_CI_groups),])
SH_CI_groups<-SH_CI_groups[order(SH_CI_groups$Year,SH_CI_groups$ITS2_14, SH_CI_groups$Community),]
SH_CI_groups #write.csv(SH_groups,"SH_Communities_multicomp.csv")
SH_Comties_I <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw+
facet_grid(~ITS2_14, labeller = labeller(ITS2_14=c(
`C_a` = "C_a (2014)",
`C_b` = "C_b (2014)",
`D1` = "D1 (2014)"))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.2),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
SH_Comties_I
#ggsave(file="Figure_3.svg", plot=SH_Comties, width=6, height=4)8b. Changes in Cladocopium abundance
CH_Com <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(Community))) + ggthe_bw+ ggtitle("a.")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.8),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.7),
name="\n Cladocopium to host cell ratio (C/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" ", " ", " "))
#CH_Com
CH_Com2ITS2<-CH_Com + facet_grid(~ITS2_14)
# CH_Com2ITS2
CH_I<-lmer(logC.SH ~ Year* Community* ITS2_14 + (1|Colony), data=data.Uva)
summary(CH_I)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * Community * ITS2_14 + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 107.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.72727 -0.28740 -0.04472 0.44178 1.54498
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.1101 0.3318
## Residual 0.3945 0.6281
## Number of obs: 55, groups: Colony, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.66014 0.29000 39.38253 -5.725 1.21e-06 ***
## Year15 -0.56741 0.36265 23.92752 -1.565 0.1308
## Year16 -0.77804 0.36265 23.92752 -2.145 0.0423 *
## CommunityMixed -0.38598 0.48624 40.68407 -0.794 0.4319
## ITS2_14C_b 0.23717 0.49523 42.78435 0.479 0.6344
## ITS2_14D1 -0.91902 0.63611 42.99226 -1.445 0.1558
## Year15:CommunityMixed 0.31035 0.70821 33.47664 0.438 0.6640
## Year16:CommunityMixed -0.34608 0.85164 32.26324 -0.406 0.6872
## Year15:ITS2_14C_b 0.55621 0.57741 27.86114 0.963 0.3437
## Year16:ITS2_14C_b 0.36653 0.83852 29.64957 0.437 0.6652
## Year15:ITS2_14D1 0.04458 0.89702 30.76732 0.050 0.9607
## Year16:ITS2_14D1 1.23468 1.06989 31.29787 1.154 0.2572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmntM ITS2_14C ITS2_14D Y15:CM Y16:CM
## Year15 -0.625
## Year16 -0.625 0.500
## CommuntyMxd 0.000 0.000 0.000
## ITS2_14C_b -0.586 0.366 0.366 -0.655
## ITS2_14D1 -0.456 0.285 0.285 -0.764 0.767
## Yr15:CmmntM 0.000 0.000 0.000 -0.628 0.411 0.480
## Yr16:CmmntM 0.000 0.000 0.000 -0.501 0.328 0.383 0.330
## Y15:ITS2_14C 0.393 -0.628 -0.314 0.546 -0.740 -0.596 -0.508 -0.278
## Y16:ITS2_14C 0.270 -0.216 -0.432 0.315 -0.470 -0.364 -0.219 -0.810
## Y15:ITS2_14D 0.253 -0.404 -0.202 0.496 -0.473 -0.610 -0.790 -0.261
## Y16:ITS2_14D 0.212 -0.169 -0.339 0.399 -0.385 -0.498 -0.263 -0.796
## Y15:ITS2_14C Y16:ITS2_14C Y15:ITS2_14D
## Year15
## Year16
## CommuntyMxd
## ITS2_14C_b
## ITS2_14D1
## Yr15:CmmntM
## Yr16:CmmntM
## Y15:ITS2_14C
## Y16:ITS2_14C 0.404
## Y15:ITS2_14D 0.655 0.260
## Y16:ITS2_14D 0.327 0.791 0.352
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients #step(CH_I)
CH_I.emmc<-emmeans(CH_I, ~Year*Community* ITS2_14)
CH_I_groups<-cld(CH_I.emmc)
CH_I_groups<-as.data.frame(CH_I_groups[complete.cases(CH_I_groups),])
CH_I_groups<-CH_I_groups[order(CH_I_groups$Year,CH_I_groups$ITS2_14, CH_I_groups$Community),]
CH_I_groups #write.csv(CH_I_groups,"CH_I_groups.csv")9b. Changes in Durusdinium abundance
DH_ComI <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(Community))) +
ggthe_bw+ ggtitle("b.")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.2),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-6, -0.7),
name="\n Durusdinium to host cell ratio (D/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
DH_Com2ITS2<-DH_ComI + facet_grid(~ITS2_14)
DH_Com2ITS2
DH_LM_ITS2<-lmer(logD.SH ~ Year* Community* ITS2_14 + (1|Colony), data=data.Uva)
summary(DH_LM_ITS2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * Community * ITS2_14 + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 66.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0953 -0.5756 0.0079 0.3844 3.6310
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.0000 0.0000
## Residual 0.1687 0.4107
## Number of obs: 58, groups: Colony, 22
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -4.86573 0.16766 48.00000 -29.021
## Year15 0.07049 0.33532 48.00000 0.210
## Year16 3.21459 0.23711 48.00000 13.557
## CommunityDurusdinium-only 0.27806 0.50992 48.00000 0.545
## ITS2_14D1 3.15614 0.23711 48.00000 13.311
## Year15:CommunityDurusdinium-only -0.35506 0.32897 48.00000 -1.079
## Year16:CommunityDurusdinium-only -0.22094 0.38416 48.00000 -0.575
## Year15:ITS2_14D1 0.36688 0.42745 48.00000 0.858
## Year16:ITS2_14D1 -3.09697 0.41068 48.00000 -7.541
## CommunityDurusdinium-only:ITS2_14D1 -0.08557 0.45916 48.00000 -0.186
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.834
## Year16 < 2e-16 ***
## CommunityDurusdinium-only 0.588
## ITS2_14D1 < 2e-16 ***
## Year15:CommunityDurusdinium-only 0.286
## Year16:CommunityDurusdinium-only 0.568
## Year15:ITS2_14D1 0.395
## Year16:ITS2_14D1 1.09e-09 ***
## CommunityDurusdinium-only:ITS2_14D1 0.853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmnD- ITS2_1 Y15:CD Y16:CD Y15:IT Y16:IT
## Year15 -0.500
## Year16 -0.707 0.354
## CmmntyDrsd- 0.000 0.000 -0.232
## ITS2_14D1 -0.707 0.354 0.500 -0.232
## Yr15:CmmnD- 0.000 0.000 0.000 -0.293 0.360
## Yr16:CmmnD- 0.000 0.000 0.000 -0.753 0.309 0.389
## Y15:ITS2_14 0.392 -0.784 -0.277 0.129 -0.555 -0.500 -0.171
## Y16:ITS2_14 0.408 -0.204 -0.577 0.671 -0.577 -0.208 -0.713 0.320
## CD-:ITS2_14 0.000 0.000 0.258 -0.900 0.000 0.000 0.558 0.000 -0.596
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## convergence code: 0
## boundary (singular) fit: see ?isSingular #step(DH_LM_ITS2)
library(emmeans)
DH_I.emmc<-emmeans(DH_LM_ITS2, ~Year*Community* ITS2_14)
DH_groups_I<-cld(DH_I.emmc)
DH_groups_I<-as.data.frame(DH_groups_I[complete.cases(DH_groups_I),])
DH_groups_I<-DH_groups_I[order(DH_groups_I$Year,DH_groups_I$ITS2_14, DH_groups_I$Community),]
DH_groups_I #write.csv(DH_groups_I,"DH_I_groups.csv")Figure 4: qPCR
Figure_4b<-grid.arrange(CH_Com2ITS2, DH_Com2ITS2,nrow=2)
# ggsave(file="Outputs/Figure_4b.svg", plot=Figure_4, width=7, height=4)10. Supplmenmtary C/H D/H dynamics
Select data
log_Log.data<-data.Uva[,c("Colony","Year.2", "C.SH", "D.SH",
"D14Community", "Community","ITS2_14")]
log_Log.data<-log_Log.data %>%
filter(Year.2 != "2015")
log_Log.data$logC<-log10(log_Log.data$C.SH)
log_Log.data$logD<-log10(log_Log.data$D.SH)
log_Log.data$logC[which(log_Log.data$C.SH==0)] <- -6
log_Log.data$logD[which(log_Log.data$D.SH==0)] <- -6Plot Black and White
Community_Changes2 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony)) +
geom_abline(intercept=0, slope=1, linetype=1, color="black")+
geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
geom_abline(intercept=1, slope=1, linetype=2, color="black")+
geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
ggthe_bw +
scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
geom_point(alpha=0.5, size=4) +
geom_path(size = 1,arrow =
arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
"10,0000C : D", "1,000C : 1D", "100C : 1D","10C : 1D", "C : D"), angle = 45)+
annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
"1C : 10,000D", "1C : 1,000D", "1C : 100D", "1C : 10D"), angle = 45)
Community_Changes2
#ggsave(file="FigAndrewsFavoriteGraph_BW.svg", plot=Community_Changes2, width=5, height=5)Community_Changes3 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=D14Community)) +
geom_abline(intercept=0, slope=1, linetype=1, color="black")+
geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
geom_abline(intercept=1, slope=1, linetype=2, color="black")+
geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts (2014)",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
ggthe_bw + theme(legend.position="bottom")+
scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
geom_point(alpha=0.5, size=4) +
geom_path(size = 1,arrow =
arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
"10,0000C : D", "1,000C : 1D", "100C : 1D","10C : 1D", "C : D"), angle = 45)+
annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
"1C : 10,000D", "1C : 1,000D", "1C : 100D", "1C : 10D"), angle = 45)
Community_Changes3
# ggsave(file="FigAndrewsFavoriteGraph_qPCR.svg", plot=Community_Changes3, width=6, height=6)Community_Changes4 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=ITS2_14)) +
geom_abline(intercept=0, slope=1, linetype=1, color="black")+
geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
geom_abline(intercept=1, slope=1, linetype=2, color="black")+
geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Main ITS2 types (2014)",
labels=c("C_a", "C_b", "D1"))+
ggthe_bw + theme(legend.position="bottom",
legend.box = "horizontal")+
scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
geom_point(alpha=0.5, size=4) +
geom_path(size = 1,arrow =
arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
"10,0000C : D", "1,000C : 1D", "100C : 1D","10C : 1D", "C : D"), angle = 45)+
annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
"1C : 10,000D", "1C : 1,000D", "1C : 100D", "1C : 10D"), angle = 45)
Community_Changes4
#ggsave(file="FigAndrewsFavoriteGraph_ITS2.svg", plot=Community_Changes4, width=6, height=6)